
# ALFF and fALFF
alff <- function(x, lf, hf, fractional=F, tr=NULL, ...) {
    if (is.null(tr))
        tr <- hdr(x)$pixdim[4]
    
    xamp <- power.niftimat(x, tr, ...)
    
    # get corresponding frequency points
    lt <- round(lf * ncol(x) * tr)
    ht <- round(hf * ncol(x) * tr)
    
    n <- length(lf)
    alff <- as.matrix(sapply(1:n, function(i) apply(xamp[,lt[i]:ht[i]], 1, mean)))
    colnames(alff) <- apply(cbind(lf,hf), 1, paste, collapse=':')
    
    if (fractional) {
        tmp <- apply(xamp, 1, mean)
        falff <- alff / matrix(rep(tmp, ncol(alff)), length(tmp), ncol(alff))
        return(falff)
    } else {
        return(alff)
    }
}

power.niftimat <- function(x, tr, ...) {
    if (class(x) != "nifti_matrix")
        stop("input must be type of nifti_matrix")

    freqs <- spec.pgram(ts(x[1,], start=0, deltat=tr), plot=F, ...)$freqs
    
    xamp <- apply(x, 1, function(y) spec.pgram(ts(y, start=0, deltat=tr), plot=F, ...)$spec)
    xamp <- sqrt(xamp)
    colnames(xamp) <- freqs
    
    return(xamp)
}

# homotopic (possibly need header knowledge)
setGeneric('homotopic', function(x, ...) standardGeneric('homotopic'))

setMethod('homotopic',
    signature(x='nifti_matrix'),
    function(x, ...) {
        # Get homotopic locations
        locs <- .homotopic.coords(x)$homotopic

        # Correlate homotopic time-series
        xcors <- apply(locs, 1, function(i) cor(x[i$lh,] x[i$rh,], ...))

        # Convert to 3D nifti space?

        return(xcors)
    })
)

.homotopic.coords <- function(x) {
    space <- tolower(sub("NIFTI.XFORM.", "", hdr(x)$qform.code))
    if (space != "talairach" || space != "mni_152" || space != "mni_305" || space != "colin")
        stop("image must be in standard space and not ", space)
    
    hcoords <- coords(x)
    
    # Filter coords to only those that have values
    # and add the vector location of these coords as a 4th column
    which.hcoords <- which(hcoords[,4]==T)
    hcoords <- cbind(matrix.ijk2xyz(hcoords[which.hcoords,1:3]), which.hcoords)
    
    # Remove any coords that are x=0
    tmp <- hcoords[,1] != 0
    hcoords <- hcoords[tmp,]
    
    # Split coords based on hemisphere
    tmp <- hcoords[,1]>0
    lh <- hcoords[!tmp,]
    lh[,1] <- lh[,1] * -1 # NOTE: making x positive
    rh <- hcoords[tmp,]
    
    # Sort
    lh <- sort.coords(lh)
    rh <- sort.coords(rh)
    
    # Filter out any coords that don't have homotopic pairs
    tmp <- apply(lh == rh, 1, all)
    lh <- lh[tmp,]
    rh <- rh[tmp,]
    
    # Fix lh (make x negative again)
    lh[,1] <- lh[,1] * -1
    
    colnames(lh) <- c("x", "z", "z", "index")
    colnames(rh) <- c("x", "z", "z", "index")
    
    return(list(lh=lh, rh=rh, homotopic=data.frame(lh=lh[,4], rh=rh[,4])))
}
